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ABSTRACT 

The characteristics of the geoid surface are quantitatively de- 
scribed; a procedure for calculating the satellite altitude is developed; 
error sources are described quantitatively, mathematically modeled, 
and evaluated in computer simulation; and procedures for maximum 
likelihood processing of altimeter data for recovery of orbit and geo- 
potential information are presented for several geopotential models. 
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DETERMINATION OF THE GEOID 
FROM SATELLITE ALTIMETER DATA 

SECTION 1 
INTRODUCTION 


Development of a capability for processing satellite altimeter data has been 
initiated. This capability will lead to an improved description of the geoid and the 
geopotential, for which the satellite altimeter is uniquely suited. The surface 
illuminated by the altimeter can be thought of as a mobile satellite tracking station 
that travels beneath the satellite and makes range measurements along the local 
vertical. 

In conventional satellite geodesy, a set of fixed tracking stations observe 
the satellite. From their measurements, the satellite orbital elements, co- 
efficients of the anomalous geopotential model, and station coordinates are 
recovered in a least-squares adjustment solution. The geoid, one particular 
surface belonging to a set of equipotential surfaces, is modeled by a refer- , 
ence ellipsoid and a set of geopotential coefficients. The reference ellipsoid 
is chosen with the use of the station coordinates. A geoid determined in 
this manner is limited in accuracy and in the level of detail presented due 
to the generally sparse distribution of tracking stations and the insensitivity 
of satellite motion to the small geoid features associated with the higher or- 
der geopotential coefficients. 

Ideally, these limitations do not apply to the use of satellite altimeter data. 

The effective number of tracking stations can be as large as is necessary, and their 
distribution is global because they are on the satellite subtrack. In addition, the 
altimeter measurement is a more direct observation of the geoid height than are 
other tracking measurements. For example, if the orbit is considered known, 
then the geoid over ocean regions may be determined directly from these meas- 
urements (References 1,2, and 3). However, the orbit is not completely known 
because of several factors, such as incomplete and inaccurate force models and 
tracking system errors. Also, the altimeter measurement is subject to error 
due both to the instrument and to dynamic sea- surface elevation effects. Orbit 
error should not degrade the determination of the relative shape of small-scale 
undulations over fairly long fractional orbit arcs, due to the nearly negligible 
effect of small-scale geopotential anomalies on the satellite motion (References 
2 and 3). However, small-scale anomalies can be expected to influence the or- 
bit over long time periods because of possible resonant conditions between the 
orbit period and the earth rotation period. 
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In theory (References 3, 4, 5 and 6), more refined determinations of the 
geoid can be made by using altimeter data to adjust the geopotential model. 

This is possible because the mathematical formulations of many geopotential 
models include parameters that can be related to the spacecraft altitude. 

Among the geopotential models that were considered are the spherical harmonic 
expansion (SHE) model, local function models (both orthogonal and non- 
orthogonal), and models utilizing area- mean free-air gravity anomalies. 

The SHE model, is convenient to use because of its universal usage and the 
orthogonality of its terms. However, a global solution on the fine grid of meas- 
urements expected from the altimeter would require the simultaneous solution 
of a prohibitive number of coefficients. 

This requirement for simultaneous solution may be relaxed by using local- 
function models, generally expressed as mathematical series, each term of 
which is nonzero only in a particular geographic area. This results in the re- 
quirement that only a small subset of the coefficients of a model need be solved 
simultaneously, constrained only by the existence of altimeter data in the geo- 
graphic domain of the local functions in question. One local function representa- 
tion (Reference 4) exhibits its local properties only in difference equation form. 
The terms of this representation are an orthogonal set and are easily related to 
spherical harmonics. 

Models utilizing surface area-mean gravity anomalies (References 5 and 6) 
are local with respect to partial sums of the total series expansion. These mod- 
els are easily converted to SHE formulations and can make use of surface gravity 
data as well as altimeter data. 

An improved geopotential description using any of the above models can lead 
to a better defined geoid in a bootstrap process by improving knowledge of the 
orbit (Reference 3). This in turn leads to better geoid determination and a better 
geopotential model. Altimeter data may also be used directly to help determine 
the orbit (Reference 3). 

At the present state of development of altimetry data processing, the im- 
portant questions to be considered concern the nature and effect of error sources 
in the altitude measurement. Although some work (Reference 7) has been done 
in this area, the questions of the amenability of these errors to filtering techni- 
ques, and the requirements to be imposed on the altimeter and its mode of oper- 
ation in order to realize an improved geoid and geopotential description remain 
•largely unanswered. 
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The results and conclusions reached to date are summarized in Section 2. 
Section 3 documents the properties of the surface to be determined, i. e. , the 
geoid. This surface is first described in general; then specific geoid features 
that are important for the design of the altimeter and altimeter data processing 
are discussed. The mechanics of the altimeter measurement are presented m 
Section 4. The influence of the altimeter and spacecraft design on the measure- 
ment are described, followed by a mathematical development of the altitude meas 
urement model for simulation purposes. This serves as a background for a dis- 
cussion of the errors affecting the altimeter measurement in Section 5. Three 
categories of errors are listed, along with applicable mathematical models. Re- 
cent work involving simulation evaluation of the effect of some of these error 
sources on the derived geoid surface is also presented. A tentative procedure 
for differential correction of the geoid, geopotential, and orbit parameters is 
outlined in Appendix A and the specific form of the equations in Appendix A for 
different geopotential models is presented in Appendix B. 


Because of the many symbols used in Sections 3 through 5, a glossary is 
provided. 
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SECTION 2 


SUMMARY 


Some of the first results obtained in connection with developing the capa- 
bility for processing data from a satellite borne radar altimeter were the de- 
scription of the characteristics of the geoid surface such as undulation amplitude, 
scale length, and slopes. These results, valuable in themselves for their impli- 
cations for altimeter design specifications, enabled the development of a more 
accurate altitude calculation- algorithm, that is, one that considered geoid slopes. 
Altimeter error sources were also identified and categorized. Mathematical 
models that were developed to describe the error sources and geopotential 
models were useful in planning the differential correction process described in 
Appendices A and B. 

Evaluation of the effects of the error sources, as represented by the study 
of simulated C-band radar tracking data in subsection 5.4, will continue. This 
will identify those error sources that are important and those that are not and 
may lead to modifications of the differential correction procedure reported in 
the appendices. Computer program algorithm based on the mathematical models 
described in this report are being developed and tested. 
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SECTION 3 


GEOED PROPERTIES 

The geoid is an imaginary surface, a particular equi potential surface of the 
geopotential field. This surface is nearly spherical, with a slight flattening at 
the poles and a pear-shaped bulge in the southern hemisphere. A finer degree 
of detail can be deduced from the magnitudes of the higher order terms in a SHE 
model of the geopotential; this reveals a nonuniform distribution of surface rough- 
ness with surface relief on the order of one part in ten thousand. The geoid is 
approximated closely by a real surface, the sea surface, in ocean areas. The 
sea surface would conform exactly to the geoid in the absence of currents, tides, 
and weather- related phenomena (Reference 3). The exact shape of this mean- 
sea-level (MSL) or marine geoid surface is dictated by the spatial variations 
in the gravitational potential of the rotating earth. Its shape is important 
for this study because ideally it is the sea surface that will be defined by the al- 
timeter measurements. 

This section examines the properties of the MSL surface, represented by the 
geoid undulation. Geoid undulation is the deviation of the geoid from the best fit- 
ting ellipsoid as measured along the ellipsoid normal. Other properties consid- 
ered include the horizontal scale length of these deviations, sea slopes, and the 
importance of these properties for the altimeter measurement. 


3.1 Undulation Amplitude 

Because the magnitude of the geoid surface relief is very small relative to 
the earth’s radius, a reference surface that closely approximates the geoid is 
customarily utilized to avoid working with small differences in large numbers. 
This reference surface, which should be easy to describe mathematically, is 
usually chosen as the rotationally symmetrical ellipsoid that best fits the geoid, 
in a least-squares sense, to the accuracy that the geoid is presently known. This 
surface plays another role in that it represents an equipotential surface of the 
"normal” potential of the geopotential field. The additional factor (the "disturb- 
ing" potential) needed to fully describe the geopotential is thus made small (Ref- 
erence 9) . 

The geopotential may be described in geocentric earth-fixed coordinates as 

W(X, Y, Z) = U(X, Y, Z) + T(X, Y, Z) (3-1) 

where W = the geopotential 

U = the normal potential 

T = the disturbing potential at the point (X, Y, Z). 
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To evaluate the magnitude of geoid undulations, one compares the geoid (de- 
scribed by a particular equipotential surface) 

W(X, Y, Z) - W 0 (3-2) 


with the reference ellipsoid 


U(X, Y, Z) = U 0 


of the same potential 


U 0 = W 0 (3-3) 

Following Heiskanen and Moritz (Reference 9) and referring to Figure 3-1, 
a point P on the geoid is projected along the ellipsoid normal onto point Q on the 
ellipsoid. This distance PQ, denoted by N, is called the geoidal undulation at P. 
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One may approximate the value of the normal potential U at P by a linear 
relation 

U P = U Q + (|t) q N $ U Q - V (3-4a) 

Since by definition 

= U + T (3 -4b) 

P P P 1 ' 


and since 



(3 -4c) 


we find by substituting Equations (3 -4c) and (3 -4b) into (3 -4a) that the geoidal un- 
dulation is given by 



Figure 3-2 (obtained from Reference 10) is a contour map ofgeoid undulation 
over the globe as derived from the SAO-M1 geopotential model. This SHE-type 
model was developed by analyzing satellite tracking data. Note that this geoid 
consists largely of long-wavelength undulations spanning 20 degrees or more 
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Figure 3-2. Level Curves of Geoid Undulation at 10-Meter Intervals - SAO-M1 Model 
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in surface distance. This is due to the exclusion of tesseral terms of de- 
gree and order larger than 8 in this SHE . * The amplitude of these relatively 
low-degree undulations is significant, however. The 169-meter (m) range 
between undulation extremes will be easily sensed by an altimeter of 5-m 
accuracy. The altitude profiles that might be seen by the altimeter over 
the SAO-M1 geoid are shown in Figure 3-3. These were generated by over- 
laying Figure 3-2 with a subtrack of the GEOS-C satellite orbit at 60 degrees in- 
clination for various values of right ascension of the ascending node. A more 
realistic geoid, that is, one which depicts shorter wavelength undulations, is 
obtained from the Detailed Geoid Computation (DGC) map (Reference 11) of the 
Puerto Rico trench area, Figure 3-4. This map was generated from 1-degree 
area mean gravity anomalies , but height contours have been hand-drawn to a 
finer scale yielding 6 arc minutes resolution near the trench. Figure 3-5 is a 
geoid profile derived from Figure 3-4 along the 67th West meridian. Note the 
sharp change in geoid undulation near 19 degrees North latitude. The 10- to 
15-m deviation between the DGC profile and the satellite geoid, drawn here for 
comparison purposes, shows the difference between a SHE of 8th degree and or- 
der for the geopotential and a geopotential model developed from 1° x 1° area 
means and represents the improvement in geoid determination that can be realized 
by a direct geoid measurement technique such as satellite altimetry. 


3.2 Spectra of Geoid Undulations 

A qualitative review of Figures 3-2 through 3-5 reveals that the largest of 
geoid undulation amplitudes are associated with long-wavelength undulations. 

This is confirmed by Figure 3-6, which shows the power spectrum of the Puerto 
Rico trench profile of Figure 3-5. Note the slight peak near the 1 degree wave- 
length on the Puerto Rico trench spectrum. This is probably due to the trench 
itself because it is about 1 degree wide in this profile. The overall shape of this 
curve shows that geoid undulation amplitude decreases rapidly as the wavelength 
of the feature decreases. If this trend to smaller undulation amplitudes is ap- 
plicable to the rest of the geoid and if it continues to higher frequencies, then it 
is unlikely that features (except unusual ones, such as the Puerto Rico trench) 
of 1 -degree wavelength or less will have undulation amplitudes exceeding 2 m. 

The shape of geoid spectra have design implications for an altimeter of 
1-m precision. To detect geoid undulations whose amplitude is near the 1-m 
precision limit, the measurement rate of the altimeter according to sampling 
theory (Reference 12), must be such that at least two measurements per 


*More recent SHE models such as the latest GSFC model (Ref. 1 1) relieve this limitation somewhat in that 
terms of degree and order 16 are used. 
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Figure 3-4. Detailed Gravimetric Geoid of the Caribbean 
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Figure 3-5. Geoid Profile, 67 Degrees West Longitude 


undulation wavelength are performed. The current preliminary design sample 
rate for the GEOS-C altimeter in one operating mode is two per second, which, 
at the GEOS-C ground speed of nearly 8 kilometers per second (kps), implies 
that undulation wavelengths smaller than 8 kilometers (km) will not be detected. 
The adequacy of this sample rate for detecting all undulations of meter amplitude 
will be reexamined when detailed geoid maps that depict the 8-km undulation 
wavelength are available. 
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The detectability of all undulations of one or more meters in amplitude im- 
poses a constraint on the design of the altimeter pulse-width, beam pattern, and 
return pulse detection process, since all these factors influence the size of the 
illuminated surface spot. Within the effective illuminated spot, all undulation 
details are averaged and no information is recovered for wavelengths smaller 
than the spot diameter. If, as has been proposed (Reference 13), some form of 
return pulse leading edge detection is used, and the altimeter system is pulse - 
width limited (see subsection 4.1), then the radius, R s , of the effective illumi- 
nated spot is given in Reference 13 as 

R s = V ctH (km) (3-6) 


where 

c = the speed of light in kilometers per second 
r = the transmitted radar pulse-width in seconds 
H = the satellite altitude (nominal) in kilometers 

For H = 1100km and r = 10~ 7 sec, 


R s = 5.74km. 

The effective spot diameter is 11.48 km and the effective spot size becomes 
the limiting factor influencing altimeter horizontal resolution. 


3.3 Sea-Surface Slopes 

Another geoid property of interest is the local slope of the instantaneous sea 
surface with respect to the ellipsoid (reference surface). This slope is related 
to the geoid undulation amplitude and scale length mentioned in subsections 3. 1 
and 3. 2, and needs explicit mention because of its effect upon the altitude meas- 
urement. This natural phenomenon will be described here and treated mathe- 
matically in Section 4. 

The sea surface can assume a slope relative to the reference ellipsoid due 
to two causes: the geoid (or MSL) slope relative to the ellipsoid (see Figure 3-3) 
and the departure of the instantaneous sea surface from MSL. A mention of the 
extremes in sea-surface slopes is sufficient to point out the problems for 
altimetry. 

Slopes due to the disturbance potential are evident in Figure 3-3, which de- 
picts geoid profiles for the SAO-Ml geopotential model. This model, which 
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describes the long -wavelength geoid undulations, exhibits slopes that range from 
0 to 5 m/arc degree, or 4. 5 x 10 -2 m/km. 

On a finer scale, much higher surface slopes can be found in localized areas, 
such as near the Puerto Rico trench (Figure 3-4). In Figure 3-4, the geoid falls 
by 11m over 1 arc degree near 60 degrees West, 18 degrees North. This re- 
sults in a slope of 9. 9 x 10 " 2 m/km. 

Dynamic processes in the ocean, such as tides, barotropic weather effects 
and long swell waves, generate slopes generally of the same magnitude as those 
due to the disturbing potential (Reference 3). On the other hand, storm surges, 
tsunamis, and ocean currents generate a sharp slope at their boundaries. Of these, 
currents are of more immediate interest because they are of a less transient na- 
ture. For example, the Gulf Stream (Figure 3-7) develops a slope on its south- 
eastern edge, due to coriolis acceleration, of 2. 26m/arc degree or 2x 10 -2 m/km. 
Similar slopes can be expected along the edges of other great currents , such as 
Kuro Shio current. 



Figure 3-7. Sea Surface Departure from the Geoid Across the Gulf Stream Between 36.5 Degrees North and 
38.5 Degrees North Near 68.5 Degrees West (Northern Section) (Reference 3) 
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SECTION 4 


THE ALTIMETER MEASUREMENT 

This section describes generally how a satellite altimeter operates, with 
emphasis on those aspects of the hardware design and operation which constrain 
the measurement of altitude. A derivation of the calculation of altitude is pre- 
sented, including a first-order correction to the sea slope. 


4. 1 Hardware Considerations 

The altimeter considered in this section is mounted on a satellite which is 
gravity- gradient stabilized such that its nominal z-axis lies within a degree or 
two of the geopotential gradient (local vertical) at the satellite. (See Reference 
14 for details of this stabilization mechanism. ) Such stabilization allows the use 
of a directional antenna for the radar altimeter. Figure 4.1, an exaggerated 
view of the altimeter measurement geometry, shows the spherical shape of the 
emitted radar wavefront. Because the transmitted beamwidth is larger than the 
nominal spacecraft libration about the local vertical, the point of first returns 
from the sea surface will lie on the normal to the sea surface through the satel- 
lite, that is, the satellite subpoint. The effective size of the illuminated spot on 
the surface is determined by the transmitted pulsewidth, the beamwidth, and the 
type of return pulse detection utilized. A pulsewidth-limited system with leading 
edge detection will be assumed; the implications of this being that the altitude 
measurement represents the weighted mean altitude over the surface area illumin- 
ated after the leading edge of the pulse intercepts the surface and before the trail- 
ing edge does so. This area is significantly smaller than that encompassed by 
the antenna beam (see subsection 3.2). As mentioned earlier, this imposes a 
limit of 11.48km on the horizontal resolution of geoid undulations, for typical 
values of altitude and pulsewidth. 

Another factor that limits the horizontal resolution of the altimeter is the 
data rate of the instrument. Sampling theory (Reference 12) requires at least 
two measurements for the shortest wavelength to be detected. The satellite will 
have a ground speed of nearly 8 km/sec. A data rate of two measurements per sec- 
ond would limit the spectral resolution to 8-km undulation wavelengths, but this 
limit is superseded by that imposed by the effective spot size. 

The configuration presented by the satellite antenna offers some protection 
against the libration of the satellite. The antenna beamwidth half angle is 1.3° 
which subtends a spot of 45. 2 km diameter at the surface. As long as the satel- 
lite effective illuminated area lies inside the beamwidth spot, i.e. the satellite 
z-axis intersects the sea surface at a point at least 5.74km inside the antenna 
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Figure 4-1. Altimeter Measurement Geometry 
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beamwidth cone at the surface, the altimeter measurement will represent the 
shortest distance between the satellite and the sea surface. When satellite libra- 
tion exceeds this limit, the altimeter measurement no longer represents the true 
altitude and sophisticated error-correcting models and procedures must be 
employed. 


4.2 Measurement Geometry 

Because the altimeter measurement is nominally the shortest distance be- 
tween the satellite and the sea surface, the measurement is along the normal to 
the sea surface that passes through the satellite. Although this direction is not 
generally parallel to the ellipsoid normal at that point (see subsection 3 . 3 and 
Figure 4-1), it is close enough for a zero-order approximation to assume that 
the ellipsoid normal is the measurement direction. At the satellite, the z-axis 
is controlled to the local vertical. Because equipotential surfaces at satellite 
altitude are even closer approximations to an ellipsoid, this local vertical can 
also be assumed to be the normal to an ellipsoidal surface. However, this nor- 
mal is not generally the same as the reference ellipsoid normal due to curvature 
of the gradient of the normal potential. 

Assuming that the measurement is made along the ellipsoid normal at S', 
the factors that influence the error in the measured altitude are represented in 
Figure 4-1. The geoid differs from the ellipsoid both in distance and in slope, 
and, to a lesser extent, the sea surface itself differs from the geoid in the same 
manner. The geoid undulation, N, is the unknown parameter to be solved for, 
but its determination is made difficult by the slope of the sea surface relative to 
the ellipsoid, the deviation of the sea surface relative to the geoid, the less-than- 
perfect knowledge of the geoid, and the uncertainty in the satellite position (orbit 
uncertainty) . 


4.3 Calculated Altitude 


To simulate the altitude measurement, it is necessary to calculate a number 
that in theory represents the actual altitude. This capability is also needed when 
forming altitude residuals, that is, the difference between observed and calcu- 
lated altitudes, for least-squares differential correction of geoid and geopotential 
parameters. The calculation of altitude is accomplished herein by an iterative 
procedure. 

Figure 4-1 shows the antenna of a spacecraft-borne altimeter sensor system 
that is oriented in the direction of the gravity acceleration vector (g) at the satel- 
lite, assuming a gravity-gradient stabilized sensor platform. The height of the 
satellite is measured in the direction of shortest signal echo, which is along the 
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normal to the instantaneous sea surface. This direction is very close to the 
local normal, n, to the geoid (i.e. the direction of line segment SP), especially 
over ocean areas. The separation of geoid and ellipsoid at P and Q, respectively, 
along the ellipsoidal normal, n', is defined as the geoid undulation, N. 

Ideally, the satellite altitude at P should be calculated along the line from 
P to the satellite, of distance H. However, the location of point P is unknown, 
so approximate techniques are used. Although N does not coincide with the SP 
line segment, the height of the satellite above the ellipsoid may be approximated 
initially as 

h - H + N (4-1) 

where H = the satellite height above the geoid at P. The coordinates of point Q 
are also generally unknown. Thus both H and N are unknown. We must have an 
approximation to h in terms of known parameters . This approximation is devel- 
oped in the following sections. 

The position of the spacecraft with respect to the geocentric coordinate sys- 
tem (Figure 4-2) is determined at any given time from the dynamics of the sys- 
tem and the estimated state vector. Because the satellite height is a scalar 
quantity, it is immaterial whether the coordinate system is inertial or earth 
fixed. However, for the definition of geopotential, it is assumed to be earth 
fixed and Greenwich oriented. The geocentric spherical coordinates of the 
spacecraft are as follows: 



These coordinates can be used to calculate the distance between the satellite and 
the ellipsoid at closest approach, that is, point S'. 
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Figure 4-2. Geocentric Spherical Coordinates 


The actual height of the satellite above the ellipsoid at S’ along the ellipsoidal 
normal is (see Figure 4-3) given by 


(4-3) 


where X e , Y e = satellite geocentric coordinates 
0 = the geodetic latitude of S' 

r = the normal radius of curvature of the ellipsoid at S' 

This expression for h' is the desired zero-order approximation to h. 

To implement Equation (4-3), one must express r N and 0 in terms of meas- 
ured or known parameters. According to Bomford (Reference 15), the radius of 
curvature, r N , of the ellipsoid at S' is given by 

a 

j— 

2 2l 2 (4_4) 

1 - (2f - f ) sin 0j 

where a = the semimajor axis of the ellipsoid and 

f = the ellipsoid flattening, both known parameters. 

Therefore, the implementation of Equation (4-3) depends only on the deter- 
mination of the angle 0. This angle cannot be calculated in closed form (Refer- 
ence 9). An iterative solution is possible, however, starting with the known 
satellite coordinates and assuming that lAs initially negligible. 



h'= 
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Figure 4-3. Reference Ellipse Measurement Geometry 


From Equation (4-3) we may write 


V 




(r N + h') cos 4 


(4-5) 
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From Bomford (Reference 15), 

Zg = { [l - (2f - f 2 )] r N +h'| sin <f 


(4-6) 


Forming the fraction and reducing we obtain 

Z 


tan 0 = 


Jx ^ + y 2 

e e 


1 - 


-1 


(2f- f*) 


N 


r N + h ' 


(4-7) 


Setting h' = 0 in Equation (4-7) , the first approximation to 0 is obtained 


0^ = tan 


-1 


if 


X 2 + Y 2 


■[(1 -2f + f 2 )] 


This value for 0 is substituted into Equation (4-4), where r Nj is evaluated. This 
in turn enables us to approximate h' using Equation (4-3). 


h ; = 


4 


x. + Y g 


COS <f 


- T 


1 


N 


(4-8) 


1 


A second approximation for 0, denoted 0 2 , may now be computed from Equation 
4-7, that is 

T - 1 


0 = tan 

2 


-1 


i 


— z r 

X + Y 
e e 


N 


1 - (2f - f ) 


1 


r + h' 
N i 1 


(4-9) 


and so on, until 0, and h'j converge. 

The zero-order approximation to the measured altitude is therefore (see 
Figure 4-3) 

H' - 6h' = h' - N' - 6h' 


where 5h' s = the deviation of the sea surface from the geoid at S' (see Figure 
4-1) 

N' = the geoid undulation at S'. 

The geoid undulation N' at the point S' is given by Equation (3-5) as 

N* (0‘, M = 
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where W = W (r se , 0' , X ) is the estimated potential S’ 


U 


■ se 


= U (r se , <p' , X ) is the normal potential S' 
= the geocentric radius to S' 


X = the longitude of S' 

4>' = the geocentric latitude at S' 

7 = the gradient of U at S' 

These parameters must also be evaluated in terms of known parameters. 


The gradient, 7 , is given by Heiskanen and Moritz (Reference 9) as 


y - V e ( x + f 2 sin 2 0 + f 4 sin%j 


where f-, = -f+ 4 -m+ 4 f 2 


26 , ,15 2 

-=- fm + — - m 
7 4 


f . = - y f 2 + ~ fm 
4 2 2 


m + ^ ^ -~- a (lo is the rotation rate of earth) 

2 Te 


Further, 


T e = 978. 049 gal 

00 = 0. 72921 151 x 10" 4 rad/sec 


The value of m may be obtained by an iterative process from 


and 


_ 00 a 3 2 

m t m . 

n y 2 n-1 


m = 0.00344986 
o 


for n = 1 , 2 , 


The geocentric radius r se and latitude f may be evaluated from 


r = r 
se N 


X 

1 - sin 2 0 (4 f - 6f 2 + 4f 3 - f 4 )j 2 
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and 


0' = tan 

To this point we have detailed the derivation of a zero order approximation 
to the altitude measurement by calculating the altitude to P'. This derivation 
assumes negligible surface slope between P and P'. Figure 4-3 illustrates the 
effect of geoid slope on the altitude calculation. The difference between H at P 
and H' at P' shows the need for a better approximation to the actual altitude H 
above the geoid at P. If one assumes a first-order or planar approximation to 
the geoid slope, the discrepancy between H and H' is found to be as much as 
2.24m for the geoid slope values presented in subsection 3.3. For example, 
from the geometry of Figure 4-4, 


1- (2f - f 2 ) 


tan 0 


(Reference 15) 
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|H-H'| = (H tan 5 b ) (tan e) 

where (tan e): The geoid slope in meters per meter surface distance 


If 5 b = 1.3 degrees 
H = 10 3 km 

and tan e = 9. 9 x 10~ 5 m/m then 


H - H' 


= 2. 24 m 

MAX 


To implement the first order approximation one must first estimate the co- 
ordinates of P. This requires an estimate of the coordinates of Q and an evalu- 
ation of geoid slope. 


Geoid slope is expressed by the angle difference (e ) between the geoidal and 
ellipsoidal normals, called deflection of the vertical. The components of e at 
P in the meridian and prime vertical planes are £ and rj, respectively, and are 
calculated through the partial differentials of geoid undulation as formulated by 
Vening Meinesz (Reference 9). 

^ = " R 

1 

17 R cos $ 3T7 . 


where 



1 

+ V ) 


(4-11) 


R = 


a Irrr 


Note that <t> is the geodetic latitude of Q. 
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Satellite 
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If one knows the geodetic coordinates of point Q (see Figure 4-3), that is, <l> 
and ft , then one can calculate the astronomical coordinates <h* , ft* (see Figure 
4-1) of P by means of 


<$>*=$+£ 


and 


ft* = ft + 


V 

cos 


The expressions for the deflection of the vertical depend upon the mathe- 
matical form of the geopotential models. If the geopotential at Pis modeled by 
a SHE, then the disturbing potential at P may be expressed as 


¥ E (f )" E (< 

n = 2 m =0 \ 


C cos m 
nm 


U + 


S si n m 
nm 


P (sin $') 
nm 


(4-12) 


where 


GM 

^nm » ®nm 

Pnm < sin *') 
L* 


is the gravitational constant of the earth 

are SHE coefficients 

are associated Legendre functions 

indicates that the terms with n=2, m=0, n=4, and m=0 are 
deleted from the summation (Reference 9) , and the geo- 
centric radius to P, denoted by r, is approximated by 


r = r se + N’ (4-13) 

The geoid undulation N'may be approximated using Brun's formula (Equation 
3-5) and Equation (4-12). Replacing d/94> by d/90' , one obtains upon substitu- 
tion into Equation (4-11), 


_l_ GM 
’ y R r 


T * n 

Efe)X 

n = 2 m=0 ' 


r p 

C cos m n+ S sin mfih- 
nm nm / a0 


nm (sin 0 ') 


(4-14) 


and 


I'wtsrr ^E(f)” E 

n = 2 m =0 


m|C sin mJ2 - S cos m P (sin 4 ' 
nm nm nm 


(4-34) 


where the partial with respect to 7 is ignored as being of second order 
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then 


dP 

nm 

(sin 

4 ') 

n 

d?<' 



COS 0' 

dP 

nm 

(sin 

4 ') 

sin 4' 

d^l' 
for m 

4 0. 


j 1 

COS 


dV [ si ^ lp „ m ( si ”?')- P n.l,m (sin H for m=0 ’ and 


The geocentric rectangular coordinates of the satellite i.e. X e , Y , Z , are 
assumed to be known. Therefore, 0 , 0', and X can be computed directly from 
Equations (4-6) through (4-10). The coordinates of P and Q, i.e. the angles $ 
and £2 and <h*and £2*, are yet unknown. If, as in Figure 4-4, the ellipsoid and 

geoid are assumed to be planar within the area subtended by the altimeter antenna 
beam, then 


and 


4 > = 0 


£2 = X 

The deflection of the vertical at P' is now the same as at P and the coordi- 
nates of P (<t>* , £2 * ) may be calculated with the use of Equations 4-11 through 
4-14. 


Knowledge of the coordinates of P permits better approximation to the cal- 
culated height above the ellipsoid. The direction cosines of the local vertical at 
P are 


cos 

a = 

cos 

(J)* 

cos £2* 

cos 

(3 = 

cos 


sin fi* 

cos 

ii 

sin 




(4-15) 


This is the direction of the true altitude measurement, assuming first-order 
(planar) corrections to the geoid slope. One can now solve for the first-order 
approximate height, h, above the ellipsoid. The equations of ellipsoid and the 
local vertical are 

b 2 X 2 +b 2 Y 2 +a 2 Z 2 ; = a 2 b 2 (4-16) 
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X 


(4-17) 


_ COS a / Z _ z ) + X 
cos 7 e e 

Y = cos &- (Z - Z ) + Y 
cos y e e 


(4-18) 


where b is the semi -minor axis of the reference ellipsoid. 


The simultaneous solution of Equations (4-16) through (4-18) for the z coor- 
dinate gives 

1 

-b 2 (A+B) _+ jb 4 (A+B) 2 - 4[a 2 + b 2 (Q 2 +P 2 )][b 2 (C+P-a 2 )]j 


+ b 2 (Q 2 +P 2 )] 


(4-19) 


where A =2Q(X - QZ ) 

* 6 6 

B = 2P (Y - P Z ) 
e e 

C = (QZ e - X e ) Z 

D = (PZ e - YJ 2 

Q = COS O'/ COS y 

p = cos /3/cos y 

The correct solution of Equation (4-19) has the same sign as Z e . Finally 

1 

2 2 

h = [(X e - X) 2 + (Y e - Y) 2 + (Z e - Z) J 

The calculated altitude measxirement to the sea surface is thus 

H - 6h = h - N (*,0, r) - 6h (4-20) 

s s 

where Sh s = the deviation of the sea surface from the geoid at P. 
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SECTION 5 


ERROR SOURCES 

The error sources that affect the estimate of the geocentric radius to the 
geoid at the measurement point are of three types: those that contribute to un- 
certainty in the geocentric distance to the satellite, those that cause the meas- 
ured surface to deviate from the true geoid, and those that affect the accuracy 
of the measurement itself. In the sections that follow, the dominant error sources 
(i.e. , those that contribute errors larger than about 10 centimeters (cm) in each 
category) are identified. * Most of these sources may be modeled by mathemat- 
ical formulae whose parameters are adjustable in a least-squares differential 
correction procedure (see Appendix A). A few sources may be treated by selec- 
tive filtering matched to their error spectra. 

After identification of the sources, their characteristics and mathematical 
models are discussed, followed by an evaluation of their effects. Error models 
developed herein are accurate to first order: second order effects such as earth 
curvature are neglected. Note that although GEOS-C satellite mission param- 
eters are used frequently as an example, the error model characteristics are 
general for all near-earth satellite altimeter missions. Additional error sources 
may exist for a particular mission or a particular altimeter design but the char- 
acteristics of these sources are not sufficiently well known at this time to allow 
their discussion in this report. 


5.1 Orbit Errors 

In general, the geocentric geoid radius is determined by the difference be- 
tween the geocentric satellite radius and the altimeter measurement. Any orbit 
error that causes a radial deviation from the expected trajectory will show up 
directly in the calculated geoid radius, or in the derived geoid undulation. 

Table 5-1 lists the factors that can cause more than 10 cm of error in the 
predicted satellite radius. All these sources are currently modeled and imple- 
mented in general orbit prediction programs, such as NONAME, GEOSTAR, and 
NAP. The correction magnitude needed represents the importance of the error 
source, not the uncertainty in altitude after correction of a particular error. 


*This choice of error threshold assumes a nominal altimeter measurement precision of ±1 m. 
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Table 5-1 

Order of Magnitude of Expected Altitude Corrections 


SOURCE 

TYPICAL MAGNITUDE 
CORRECTION NEEDED 

TROPOSPHERIC REFRACTION 

2 M (REFERENCE 17) 

IONOSPHERIC REFRACTIVITY AND 
REFRACTION 

10 CM 8 

TROPOSPHERIC REFRACTIVITY 

1 M 

TRACKING STATION POSITION ERROR 

10 M (REFERENCE 10) 

MEASUREMENT BIAS 

10 M (REFERENCE 7) 

TIMING BIAS 

10 M 

PHASE ANGLE ERROR (FOR PHASE 
SYNCHRONOUS TRACKING SYSTEMS) 

10 M 

MEASUREMENT SCALE ERROR 

1 M 

DYNAMIC MEASUREMENT TIME LAG 

10 CM 

ATMOSPHERIC DRAG 

1 CM b 
TO 1 KM 

SOLAR WIND 

b 

1 CM 
TO 1 KM 

GEOPOTENTIAL ERROR 

20 M (REFERENCE 7) 


8 DEPENDS ON TYPE OF MEASUREMENT. 

^DEPENDS ON ARC LENGTH, SATELLITE, AND ALTITUDE. 


These numbers depend very much on the type of tracking system used. For ex- 
ample, a laser tracking system that may suffer a 1-m error due to ray path 
bending in the troposphere (see Ref. 18, pp. 21-23) is completely unaffected by 
the ionosphere. On the other hand, Doppler tracking usually incorporates a 
built-in correction for ionospheric effects, and resulting ionospheric errors are 
on the order of centimeters. 

Tropospheric refraction is the time delay experienced by electromagnetic 
signals due to the greater-than-unity value of the tropospheric index of re- 
fraction. As mentioned earlier, refractivity refers to ray path bending. The 
amount of bending (and the consequent range error) varies with radio fre- 
quency. Tracking station position error is mislocation of the station in the 
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geocentric coordinate system. Measurement bias is merely a calibration error 
in the tracking system. Constant error in the estimate of the time of the meas- 
urement is called timing bias. Phase-angle error is analogous to a measure- 
ment bias for such systems as Minitrack that rely on measurement of phase an- 
gle. Scale error is an error whose magnitude is proportional to the magnitude 
of the measurement. Dynamic measurement time lag refers to the satellite's 
change in position between the time that the tracking signal is sent and the time 
that the measurement information is received at the station. The radial position 
error due to dynamic time lag is proportional to satellite radial acceleration. 
Atmospheric drag tends to cause a secular decrease in the value of satellite apo- 
gee (and orbit radius) leading eventually to atmosphere reentry of the satellite. 
Solar wind is a periodic force (once per satellite revolution) leading to changes 
in orbit eccentricity and hence a change in satellite orbit radius at both apogee 
and perigee. Both of these sources of radial error are included in the force 
model of the general orbit prediction programs and are included in the corres- 
ponding least-squares correction process. Error in the geopotential can cause 
both short- and long-term radial errors, although the dominant frequency in- 
volved appears to be that of the revolution of the satellite about the earth (Refer- 
ence 7). 


The model for tropospheric refraction is in the form of a correction to the 
measurement. For example, for range, p , the correction due to refraction 
(Reference 16) is 


AP = 


2r ( 
s s 


n - l) 2 [ 1 + n s 4(n -1) ] 


(meters) 


L . 2 ? 

i in E + (i + sin : E) 


(5-1) 


where r s = the geocentric radius to the station 
n s = the index of refraction 
j?2 = 28800/r s 

E = the elevation angle of the ray path from the station to the satellite. 


Other more sophisticated models for tropospheric refraction correction such 
as that by Hopfield (Reference 18) may also be used. 


Refractivity correction can be implemented by means of a correction to the 
elevation angle as follows (Reference 16), 


2 (n - 1) cos E 

A E = — = — (radians) 

sin E + (i + sin E) 


(5-2) 
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The models for ionospheric effects are similar in form to Equations (5-1) 
and (5-2), for a fixed frequency tracking signal. 

Corrections to range measurements due to other sources are modeled (Ref- 
erence 16) as 

<5p = p x B E + P y B N + p z B u + B + P T + Ap r f R + t & + p S + pL 
s y s s 

(5-3) 

where p x , p y , p z = the partial derivatives of range with respect to the com- 
ponents of tracking station error X s , Y s , Z s 

B e , B n , By = the coefficients of tracking station error in the East, 
North, and vertical (up) directions 

= the measurement bias coefficient 

= the time rate of change of range 

= the coefficient of the timing bias error model 

= the specific ionospheric refraction correction 

= the coefficient of the ionospheric refraction 

= the model for specific phase error 

= the phase-error model coefficient 

= the coefficient of the range scale -length error model 

= the coefficient of the range dynamic lag error model 

An evaluation of the effect of some of these error models is presented in 
subsection 5.4. 

In the Network Analysis Program (NAP) general orbit prediction (Reference 
16), atmospheric drag, solar wind, and geopotential effects are all evaluated in 
a differential correction process in accordance with the mathematical models 
described in Reference 16. The effect of error in the geopotential model, illus- 
trating the situation in which geopotential correction is not applied, is depicted 
by Marsh and Douglas (Reference 7) . An error in spherical harmonic coefficients 
of the difference between coefficients of two contemporary SHE geopotential mod- 
els was simulated and its effect, on a typical GEOS-C satellite orbit was plotted. 
Radial errors in orbit prediction of 20 -m amplitude with a period equal to the 
orbital revolution resulted. This indicates that when the geopotential error is 
merely modeled without differential correction, the effective error in radial or- 
bit prediction may be removed from the altimeter measurement residuals by a 
filter set at the satellite revolution frequency. 


B 

P 

T 

Ap rf 

R 
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The effect of small-scale anomalies in the geopotential, such as the Puerto 
Rico trench, are typically not modeled by SHE models. This trench causes a 
perturbation of about 200 milligals (mgal) in surface gravity acceleration over an 
area 1 degree wide and 8 degrees long. To study the effect of such anomalies 
on the radial component of the GEOS-C orbit, a concentrated mass (mascon) was 
located at a depth of 100 km and its effect on the orbit was simulated in the NO- 
NAME orbit prediction program. The magnitude of the mass was chosen to yield 
the observed surface gravity anomaly (Reference 11) . The short-term effect on 
the radial component of the orbit is negligible compared to the 20 m error due to 
uncertainty in lower degree and order terms of the geopotential SHE. The radial 
deviation of the orbit from the nominal orbit is plotted in Figure 5-1. The true 
effect of the local anomaly is evident within minutes of the time of closest ap- 
proach to the mass concentration. The long-period phenomenon, with a frequency 
equal to the satellite revolution, is erroneous and is due to the aliasing of the ef- 
fect of the- added mass with the low-order zonal terms of the SHE geopotential 

model. 


Radial Orbit 
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5.2 Hardware Errors 


It is evident that any error in the altimeter measurement directly affects the 
accuracy of the determination of the geoid radius, inasmuch as the geo id radius 
is calculated from the altimeter measurement. Table 5-2 lists the factors that 
cause altimeter measurement error in excess of 10 cm. These factors are de- 
pendent on many parameters and can cause a wide range of error values. Typical 
expected values for these factors are presented in Table 5-2. 


Table 5-2 

Altimeter Measurement Corrections 


SOURCE 

TYPICAL MAGNITUDE 
CORRECTION NEEDED 

ALTITUDE RATE 

30 CM Note 1 

ANTENNA OFFSET 

1.5 M Note 2 

SPACECRAFT LIBRATION 

600 M Note 3 

TRACKING LOOP NOISE 

75 CM (Ref. 17) 

CLOCK INSTABILITY 

10 CM (Ref. 17) 

CALIBRATION ERROR 

1.2 M Note 4 

ALTIMETER DRIFT 

30 CM Note 5 

REFRACTION 

25 CM (eqn 5-8) 

REFLECTIVITY 

60 CM (Ref. 19) 


Note: 1. From consideration of orbit ellipticity (e < 0.006), geoid 
slopes (Section 3.3) and a pulse rate of 6 per minute. 

2. From geometrical considerations (eqn. 5-5) at 1.4° offset 
angle. 

3. From geometrical considerations only (eqn. 5-5) at 3.3° 
libration angle. 

4. Assuming a 4 nsec receiver delay time. 

5. From equation 5-7 for assumed drift rate of 3cm per 100 
hours of operating time and an operating time interval of 
1000 hours. 
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Altitude rate is analogous to the dynamic tracking lag discussed in subsec- 
tion 5.1. It is due to the time rate of change of altitude of the satellite and the 
finite time lag between transmission and return of the altimeter pulse. The er- 
ror in altimeter reading due to altitude rate may be modeled as 

"h. - h. 

Ah a = h£ = ^ £ (5-4) 

i+l i 

where ii = the time rate of change of altitude rate 
h, = the altitude rate at time t; 

= the undetermined coefficient of the error model (sec 2 ) 

Antenna offset refers to the misalignment of the antenna beam axis and the 
spacecraft stabilized axis. If the antenna is offset from the stabilization axis by 
an amount larger than the effective angular baamwidth, then the radar return to 
the spacecraft will not originate from the closest surface point to the satellite, 
but from some point whose distance is greater than the altitude. The altitude 
error due to this gross offset depends on the beamwidth half angle and the offset 
angle, and may be expressed for offset less than 5° from simple geometry by 

/a 0 - 5 b \ 

Ah 0 = H tan l ) tan (a 0 - 5 b )fora 0 > 5 b . (5-5) 

where a 0 = the offset angle 

5 b = the beam half angle 
H = the nominal altitude 

Note: this formulation assumes that sufficient signal strength exists at high 
libration angles for normal operation of the altimeter signal detector. 

Even when the offset is not enough to cause the closest point on the sea sur- 
face to fall outside the beam pattern, another type of antenna offset error is sig- 
nificant. This is the unequal distribution of transmitted radar signal around the 
measurement point that occurs when the measurement point lies off of the beam 
pattern symmetry axis. The absolute magnitude of the signal strength at this off- 
axis point is also generally different from the strength on-axis, thus causing er- 
rors in leading edge detection (Reference 17). At the present time, not enough 
is known about the design details of the altimeter to warrant an explicit error 
model for this comparatively minor antenna offset effect. 
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Spacecraft libration can cause errors whose characteristics are similar to 
those described above as antenna offset errors. In this case, the satellite at- 
titude is offset from the desired stabilization direction due to orbit eccentricity, 
thermal effects, or solar wind (Reference 14). The offset is usually periodic 
in nature with a frequency equal to that of the satellite revolution about the earth. 
Replacement of the angle a 0 in Equation (4-24) by 

V A i sin <“ 0 t+ V (5 -6) 

where = the spacecraft libration angle 
Aj i - libration model coefficient 
oj q = the orbital angular rate 
t = time from satellite epoch 
<t>jl = the libration phase angle 

yields a satisfactory error model for libration effects. 

The altimeter is subject to many sources of error inherent in the electronic 
equipment (Reference 17). The exact error models developed to correct for these 
errors must await the design and development of the actual altimeter. Some of 
the more important sources are instrument delay, thermal noise, tracking loop 
jitter, clock instability, instrument delay variations, and data readout ambiguity 
(Reference 17). Instrument delay contributes a bias to the measurement and is 
included in Table 5-2 as calibration error. 

Clock instability is a random error due to short-term d-’dft instability in the 
altimeter's internal clock. The remainder of the instrument error sources con- 
tribute random errors and are lumped together as tracking loop noise in Table 
5-2. 

Altimeter drift is an error that is linear with time. Caused by such effects 
as electronic component aging, it is modeled by 


Ah = D(t - t ) 
D o 


(5-7) 


where D = the undetermined error model coefficient 

t-t = the time interval from present time t to initial time t 0 
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Atmospheric refraction also affects the altimeter measurement, although 
not as much as it does the tracking measurements. The vertical altimeter signal 
path through the atmosphere traverses less air mass and is not subject to bend- 
ing. Hopfield's recent tropospheric refraction model (Reference 18), developed 
for Doppler measurement correction, has been adopted for correction of the al- 
timeter measurement. Hopfield's model, with vertical incidence of the altimeter 
signal on the surface, reduces to 

Ah n =2x10 (43. 13 . 5. 206 sin + 12 N'j'wJ (5-8) 

where N xd and N Tw = refraction coefficients for the dry and wet components 

of the atmosphere (N xd = 265 and N Xw = 110 in mid- 
latitude marine summers). 

0 = geographic latitude 

Variations in reflectivity of the sea surface within the illuminated altimeter 
footprint can cause a distortion of the altimeter return pulse with consequent al- 
timeter error. Pierson and Mehr (Reference 19) have found that the altimeter 
bias caused by the variation of reflection statistics with wave height can amount 
to tens of centimeters. The model describing the effect of reflection statistics 
on altimeter error requires knowledge of the transmitted pulse shape and detec- 
tor design. Implementation of this model for evaluating and correcting surface 
reflectivity error must await a firm altimeter design. 


5.3 Sea-Surface Errors 

The instantaneous sea surface generally lies within a meter or two of the 
geoid (Reference 3), but such deviations cannot be neglected if an altimeter ac- 
curacy of 1 m or less is desired. The altitude measurement represents the 
shortest distance between the satellite and the sea surface. The deviation of the 
sea surface from the geoid is denoted 5h s and is defined along the ellipsoid nor- 
mal (see Figure 4-1). Phenomena that can cause 5h $ to be 10 cm or larger are 

listed in Table 5-3. 

Referring to Table 5-3, swell describes long- wavelength water waves of the 
sort that propagate thousands of miles from the storm that spawned them. They 
typically have a sinusoidal shape, awavelengthof 200mto2km, and propagate prin- 
cipally in one direction. Wind waves differ from swell in that they are shorter in 
wave-length, not obviously directional, and of irregular shape. Usually associated 
with strong winds, they have wavelengths from 1 cm to 300 m. Strong winds 
sometimes pile up water against the coast causing a local rise in the sea surface. 
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Table 5-3 


Sea Surface — Geoid Delation Sources 


SOURCE 

TYPICAL , 
MAGNITUDE 1 

SEA SWELL 

1 M 

WIND WAVES 

1 M 

STORM SURGE 

10 CM 

BAROTROPIC DEPRESSIONS 

10 CM 

CURRENTS 

1 M 

TIDES 

1 M 


''reference 3. 


This effect is called storm surge. Barotropic depressions caused by high- and 
low-pressure areas of global weather patterns result in a local sea-surface de- 
viation from the geoid. 

Ocean currents can cause local sea-surface deviations due to the coriolis 
effect. Consequently, the observed surface deflection, with a maximum at the 
maximum velocity point, is greatest in the higher latitudes and for East-West 
currents. 

Tides, so well measured near the coasts of continents, are nearly undeter- 
mined in midocean. Theories for the time-dependent sea-surface deviations due 
to tides have been developed (Reference 3) and the frequencies of the various tidal 
components are known. A localized time -dependent, tidal deflection error model 
has been developed for the component frequencies whose amplitudes are domi- 
nant (Reference 20). Accordingly, the expected tidal surface deviation at latitude 
(in degrees), longitude X (in degrees), and mean solar time (at Greenwich in 
hours), t, is 

Ahrp = 0.2667 

+ sin 2 $ |0K 1 cos (i 5. 041t+ 65 +X)I (5-9) 


(meters) 


+ cos $ M_ cos (28. 984t + 135 + 2X ) 


(1-3 sin $) A t + M f cos (1. 098t + 41.2) 
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where A t , M f , OK, , and M 2 are coefficients to be adjusted. 

This model was developed from the known dominant monthly, daily, and 
half -daily tidal components, aggregating components of common periods by 
choosing an average phase angle. 

The other sea-surface effects mentioned above may be modeled by a local 
sea-level bias. 


5.4 Simulation Results 

The effect on the predicted satellite position of bias in the satellite tracking 
measurement has been partially evaluated in a computer simulation. The NAP 
program was used to simulate biased and noisy C-band radar tracking data for a 
typical GEOS-C orbit. The NAP program simulates the tracking data and com- 
putes a least-squares differential correction to the GEOS-C state vector. This 
adjusted state vector represents the best estimate of the satellite state vector 
given the erroneous tracking data. The difference between the radial components 
of the satellite position vectors predicted by the simulated data and by the uncor- 
rupted data (represented by the nominal or a priori state vector) is evaluated and 
plotted in the NONAME Program (RV compare option) . This vector component 
difference represents the error in calculated altitude due to the biased tracking 

data. 

To date one simulation has been carried out: for a 6 -hour orbital arc with 
perturbed tracking data from eight C-band tracking sites. The a priori meas- 
urement bias was 6 m with a standard deviation of 3 m, which is the expected 
range bias uncertainty in accordance with officially documented tracking system 
performance (Reference 21). The stations tracking GEOS-C were assumed to be 

located at: 

Kauai, Hawaii 

Pt. Mugu, California 

White Sands, New Mexico 

Eglin AFB, Florida 

Cape Kennedy, Florida 

Grand Bahama Island, Bahama Islands 

Antigua Island, Lesser Antilles 

Pretoria, South Africa 

The elevation limits imposed precluded tracking below 5 degrees elevation 
for these stations. The average pass length was 15 minutes, and the 6-hour arc 
length allowed about two consecutive passes per station. The gravity field model 
used was the SAO-M1 field. 
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Simulated data were generated every 5 minutes during the visibility period 
for each station. The resulting predicted orbit deviates from the nominal (er- 
ror free) orbit and the radial deviations are plotted in Figure 5-2. 


Radial Orbit 
Deviation 
From Norminal 
(meters) 



The maximum radial deviation, or altitude error, due to tracking system 
bias in this case is about 1 m. This is the magnitude of the error in the altitude 
calculation due to bias and noise in C-band tracking for the situation tested. The 
altitude error is expected to vary as data rate, station distribution, and tracking 
duty cycle are changed. 


5-12 



APPENDIX A 


SAMPLE DIFFERENTIAL CORRECTION PROCEDURE 

The following discussion is provided as an example of the differential cor- 
rection procedure that might be used in processing altimeter data. This proce- 
dure involves the formation of the difference (residual) between the observed 
and calculated (expected) altitudes and the solution for differential corrections 
to the parameters of the problem, such that the weighted sum of the squares of 
the residuals is minimized. 

The residuals are formed using the difference between the expected and the 
observed altitude measurements. The observed altitude is the altitude meas- 
urement at time t r , denoted O r . The expected altitude to the sea surface may 
be obtained from Equation (4-20). To this must be added certain error models 
that enable one to express the expected altitude measurement. In functional 
form, the expected altitude measurement function is 


H = h (eT , t) - N(CT ) - <5h (T, S' , 0, A, t) - Ah (m , 0, t) (A-l) 

o o s o o o 

where h = the calculated altitude to the reference ellipsoid 

E 0 = the vector of orbital elements, and other arc-dependent para- 

meters, such as drag and solar pressure 

t = time 

G 0 = a vector of geopotential coefficients 

N = the geoid undulation function 

5h s = the sea-surface deviation function 

= a vector of tidal error model coefficients 
S 0 = a vector of local sea-surface biases 

0 = the geodetic latitude of the altimeter measurement point 

A = the longitude of the altimeter measurement point 

Ah = the functional representation of altimeter measurement errors 


m 0 = a vector of altimeter measurement error coefficients 


Written in this form, Equation (A-l) is not the best estimate of the true 
altitude at time t but is the best estimate of the actual altitude measurement, as 
affected by altimeter error sources. Note that for simplicity the expected 
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altitude to the ellipsoid, h, is not here expressed as a function of geopotential 
coefficients, G 0 . The measurement residual is 


(O-*) 

r 


where r denotes the r th measurement. 


Mathematically, this residual is equivalent to the difference function, 


A *r = 


h=l 



L 

dE h + J, 


5 

dG k + 2, 

dM 

dE, 

dG, 

dT 

h 

r k=l 

k 

r 1 = 1 

i 


dr. 


^ d>f 

+ ^S“ 

J = 1 J 


dS. + 


i = l 


~5 


m . 

1 


dm. + 60 
1 


(A -2) 


where L, q, and p are the total number of geopotential, sea-surface bias, and 
altimeter measurement error model coefficients, respectively. In matrix form, 
Equation (A-2) is 

(O - 3+7 = [A ] dE + [C ] dG + [F ] dT + [K ]dS + [L ] dm + T (A-3) 
o o o o o o ' 7 

-i. 

where 5 0 is the N-dimensional vector of unmodeled errors between O and M 
for the N measurements, r = 1, 2, . . ., N. 

To this observation equation we add the observation equation of the additional 
satellite tracking data available from satellite-to-satellite tracking and various 
ground-based tracking systems, such as laser, C-band, and Doppler, 

(O’ " M) = [Q] dE + T (A _ 4) 

The resulting observation equation may be further reduced in matrix notation to 


R = [m] Z + 


( A-5) 
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For a least-squares solution, we seek a vector Z such that 

[g] Tj = MIN (A -6) 

where [ G ] = the inverse variance/covariance matrix of the observations; that 
is, the weighted square of the unmodeled error, is minimized. 

These measurements are usually assumed to be uncorrelated so that Gy 
0, i ^ j, and the measurement weight matrix [ G] is diagonal. 

The minimization of t[ [ G ] 5, is accomplished by the solution of the nor- 
mal equation 

[m] t [g][m] z= [m] t [g]1r (A" 7 ) 


Equation (A- 7) is solved for the differential corrections Z, which are then 
applied to the error models coefficients and force model parameters. The ex 
pected altitude measurement * is recalculated based on these improved co- 
efficients and parameters, and a new residual vector R, is calculated. Equa k> 
(A- 7) is solved again, if necessary, and each iteration reduces the residuals R 
until they and the calculated differential correction are small enough that the 
parameter vectors in Equation (A-l) can be said to have converged. 

The direct solution of Equation (A-7) is straightforward but not always 
desirable. It may happen that the a priori parameter values used m forming 
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the matrix M are so different from their converged values that Equation (A- 7) 
will not lead to a converging solution. If this is the case, and if the normal equa- 
tion matrix [M] T [G] [M] is such that off-diagonal elements common to different 
subsets of parameters are small compared with the diagonal elements, then the 
parameter set Z, and the normal equation matrix may be divided into subsets to 
be solved separately. Thus an improved set of values for parameter subset 1 
may be obtained for use in solving for parameter subset 2, and so on, until a 
complete set of starting elements for solving Equation (A-7) is obtained. Of 
course, if the off-diagonal elements of [M] T [G] [M] are zero for a given subset 
of parameters, then that subset of parameters may be solved completely inde- 
pendently of the other subsets . 

This technique is currently used in processing satellite tracking data where- 
in orbit parameters are adjusted and converged separately and only one iteration 
thereafter is necessary to obtain good geopotential coefficient adjustments. A 
similar situation may exist in the case of altimeter data processing. The orbit 
parameters and geopotential coefficients should have small off-diagonal elements 
in common in the normal equation matrix. Similarly, sea-surface deviation 
parameters are relatively independent of the other parameters. On the other 
hand, altimeter measurement error parameters are closely related to orbit 
parameters in their effect on measurement residuals. Therefore, solution of 
dE and dm (hereafter called the orbit subset), separately from d T and dS (the 
surface subset), and from dG, should be fruitful. 

This particular partition of altimeter parameters is also desirable because 
of the different characteristics of these parameters. The orbit is best solved 
for using time-ordered sequences of data such as tracking station passes and 
altimeter turn-on periods, which are independent of geographical area. The sur- 
face subset, because of the localized nature of the sea-surface deviation phenom- 
ena, is best solved for using the measurements that fall in the localized area in 
question. 


Solving first for the orbit subset, Equations (A-3) and (A-4) may be written 
as 


(O -M) 


A | L 
o | o 

\ T o 

~' rr ~ — 


Q j 0 

dm f < A -*> 

L J L ° _ 

(O' - M) 


or 


R = [Ml 
L o 



(A- 9) 
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* 



where the subscript o indicates orbit subset. The corresponding normal equation 
solution for the orbit subset of parameter adjustments is 


(- ?-) ■ 

[M]J [G] [M] c 




-1 


[M] o [G] R 


(A-10) 


The updated orbit subset parameters are 


E, = 

E + dE 

1 

o 

m, = 

m + dm 

1 

o 


(A-ll) 

(A-12) 


At this time, the altimeter residuals may be processed through a high-pass 
filter to remove the low-frequency radial orbit perturbation contributed by low- 
order and degree geopotential terms. The altimeter residuals now contain high- 
frequency geoid undulations contaminated by sea- surface deviations and leftover 
uncorrected errors in the orbit subset. 

The relevant equations for solution of these orbit and sea-surface effects 
are as follows: 

j* = h" (E. , "G , t)-N(G )- 6h (t , S' , 0, X, t)- Ah(m , 0, t) (A-13) 

I ' 1 o osoo 1 


and 


(0-3^) = [F Q ]dT+ [K Q ]dS + [A J ]dE+ [L^dml & z (A-14) 




[m]^ [g] R 

os i 


(A-15) 


*2 

II 

ML 

+ 

dE 


™2 

= “i 

+ 

dm 

(A-16) 

T* 

1 

= r“ 
o 

+ 

dr 


r i 

II 

U1 

o 

+ 

dS 
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where the subscript os indicates the orbit and sea-surface subsets and [A t ] and 
[ L ! ] are the updated submatrices of partials for the orbit subset of parameters. 

The residuals that remain after the expected altitude measurement # is. re- 
calculated contain geoid undulations and uncorrected orbit and sea-surface effects. 
Equations (A-3) through (A-7) may be updated and summarized: 

*11 = ShgfTj.S^*. A. t)-Ah(m r 0, t) (A _ 1?) 

( °"*II ) = [C Q ]dG + [A 2 ]d*E+ [L 2 ]dm+ [F 1 ]dr+ [K 1 ]dS+T 3 (A-18) 

z = ^[m]j T [g] TmJj J [m]^ [g] E n (A-19) 

where the subscripts on matrices and vectors reflect the updating of the param- 
eters through previous subset solutions. 
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APPENDIX B 


GEOPOTENTIAL MODELS 

Several different mathematical models of the geopotential have been examined 
for use in the altimeter data processing simulation. These models and their re- 
lation to the calculated altitude measurement are presented below. 

B. 1 S pherical Harmonic Expansion 

In this model, the disturbing potential T is written as 


T(r, 0, X) 


GM 

r 


L* 


E E (t 


L n=2 m=0 


(C 


cos mX + S sin mX) P (sin 0 ) 


(B-l) 


The relation between T and the calculated altitude measurement # is expressed, 
using Brun’s formula (Equation (3-5)) as 


* = h (C , S ) - 6h - Ah 

y nm nm s 


(B-2) 


The appropriate partial derivatives in the altitude measurement condition 
equation (see Appendix A) are 


where 





as 


nm 



as 

nm 


(B-3) 
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and 


nm nm 


8X _ 8N GM /a \ m 

"as — = as = ~T~ \~) sinmXP n (sin r) 


B.2 Local (Sample) Function Model 

In this model, the disturbance potential is expressed in terms of sample 
functions that are nonzero only in certain geographic areas. 


T(r, 0,X) = GM 


G. ](. (r, <b, a) 


(B-5) 


where Jr j(r,0,X) are the sample functions (Reference 4). 
The calculated altitude measurement is written as 


X = h - — (G) - 6h - Ah 

y s 


and the conditions equation partial derivatives become 


(B-6) 


_ k _ _J_ _9T_ _ GM . 
°ki 9 G J i 7 9G i 7 * 
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where = 0 unless the k th observation lies in the l th geographic area. 

B.3 Surface Layer Model 

J 

This model presents the disturbing potential as due to a thin layer of 
material of varying density Gj deposited on the reference surface (Reference 22). 
The disturbing potential is expressed as 


T(r, 0, X) = 


i 


(B-8) 
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where Aq = the i th finite surface element 

i j = the distance between the point (r, <j> , A ) and the centroid of the i*h 
surface element 


The calculated altitude measurement is given by Equation (B-6) and the 
condition equation partial derivatives are 


C 

°ki 


J_£T _ 1 f f da 

a Gi -r 9 Gi -T JJ ^ 


(B-9) 


Vinti (Reference 23) also describes a surface layer model. 


B.4 Free-Air Gravity Anomaly Model 


This model (Reference 24 and 11) express the disturbing potential at a given 
point P (r, 0 , A) by a sum of gravity anomalies in a small area around P. The 
disturbing potential is written as 


T(r, p, X) = T^r, 0, A) 


+ 



Ag (0.* A •) S {*) Aa 

J J j 
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where T L 


Ag(0j, Aj) 


the disturbing potential calculated by a set of low-order and 
degree spherical harmonics which are not adjusted 

the mean free-air gravity anomaly referenced to the equi- 
potential surface of T x in the j th element of area 


* 

SW 


= the angular arc between (0,A) and , Xj ) 
= Stokes' function 
= the size of the elemental areas 


The calculated altitude measurement is given by 


>HP) = h (P) - -y (G) - 6h (P) - Ah 


(B-ll) 
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where 


/ Ag (0,.Xj) \ 


G = 


Ag (0..X.) 

J J 


Ag (0 L , X L ) 


and the corresponding condition equation partial derivatives are 


: o ki = t s( V a<7 * 
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Symbol 

[A 0 ] 

A je 

A t 

a 

a 

a Z 

(Lq 

B 


b 

P 

[ C °lJ 

c 

7 

& 
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GLOSSARY OF SYMBOLS 
Definition 

Matrix of orbital element partial derivatives 

Libration coefficient 

Tidal model coefficient 

Semimajor axis of the reference ellipsoid 

Angle between the local vertical at p and the x-axis of the geo- 
detic coordinate system 

Satellite libration angle 

Antenna offset angle 

Measurement bias 

Coefficient of tracking station error 

Coefficient of tracking station error 

Coefficient of tracking station error 

Semiminor axis of reference ellipse 

Angle between the local vertical at p and the y-axis of the geo- 
detic coordinate system 

Matrix of geopotential coefficient partial derivatives 
Speed of light (kilometers per second) 

Angle between the local vertical at p and the z-axis of the geo- 
detic coordinate system 

Gradient of U at Q 

Phase error 

Altimeter drift error coefficient 
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Symbol 


Definition 


Ag 

Ah a 

Ah D 

Ah 0 

Ah n 

Ah T 
A P 
Ap rf 
5b 

5h s 


5h s - 





5p 


Aa j 



E 


e 


t F ol 

f 


Mean free air gravity anomaly 

Altitude correction due to altitude rate 

Altimeter drift errbr 

Antenna offset altitude correction 

Tropospheric refraction altitude correction 

Tidal altitude correction 

Range error due to refraction 

Specific ionospheric refraction correction 

Half angle of altimeter beamwidth 

Sea- surface deviation from the geoid along reference ellipsoid 
normal at Q 

Sea-surface deviation from the geoid at S’ (acceleration of normal 
gravity) 

Vector of unmodeled errors 

Correction to range measurement due to error source 

Incremental surface element 

Vector of orbital and arc-dependent elements 

Elevation angle of ray path from tracking station to satellite 

Angular difference between geoid and ellipsoid normals (deflec- 
tion of the vertical) 

Matrix of tidal model partial derivatives 
Flattening - of the reference ellipsoid 
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Definition 


Symbol 


[G] 

GM 

G 

g 

* 

H 

H' 

H" 

h' 

h 


h l> h i 

hj 

ii 

V 

[K 0 ] 

[L 0 1 

L;L* 

X 

[M 0 J 


Inverse variance/covariance matrix 
Gravitational constant of the earth 
Vector of geopotential coefficients 

Gravity acceleration vector (local vertical) at the satellite 
Expected altitude measurement 
Nominal satellite altitude (km) 

Satellite altitude above the geoid along ft at S* 
Approximation to H at P 

Altitude of the satellite above the reference ellipsoid at S’ 

Approximation to the satellite altitude above the reference 
ellipsoid at P 

Successive approximations to h 

Time rate of change of altitude at time tj 

Time rate of change of altitude rate 

Component of the deflection of the vertical 

Matrix of local sea-surface partial derivatives 

Matrix of altimeter error model partial derivatives 

Altitude rate error coefficient 

Range dynamic lag; number of geopotential coefficients 
Longitude 

Matrix of parameter partial derivatives 
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Definition 


Measurement vector of "other" data type 
Tidal model coefficient 
Tidal model coefficient 

Vector of altimeter measurement error coefficients 
Order of spherical harmonic term 
Geoid undulation (meters) 

Geoid undulation at S’ 

Dry and wet refractivity coefficients, respectively 
Degree of spherical harmonic term 
Unit vector along the ellipsoid normal 
Index of refraction 

Component of the deflection of the vertical 

Observed altitude at time t r 

Tidal model coefficient 

Observation point on the geoid surface 

Number of altimeter measurement error coefficients 

Point on the reference ellipsoid obtained by projection from P 
along the ellipsoid normal at P 

Number of local sea- surface bias coefficients 

Vector of altimeter residuals 

Coefficient of ionospheric refraction 

Effective illuminated surface spot radius (km) 



Symbol 


Definition 



S 

sw 


S' 


T 

T i 

t 


r 

r o 

U 

$ 

<J>* 


Geocentric radius to P 

Geocentric radius of the satellite 

Normal meridional radius of curvature at S’ 

First approximation to r 

Geocentric radius to tracking station 

Geocentric radius to point S' on the reference ellipsoid 

Range from tracking station to satellite 

Time rate of change of range 

Partial derivative of range with respect to parameters p 
Range scale-length error 
Stokes' function 

Satellite subpoint on the reference ellipsoid 
Vector of local sea-surface biases 
Disturbing potential function; timing bias 
Disturbing potential function 
Time 

Transmitted radar pulsewidth (seconds); specific phase error 

Vector of tidal error model coefficients 

Normal potential function 

Geodetic latitide of point Q 

Astronomical latitude of P 
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Symbol 


Definition 


0 


0 1 > 02 5 


0 

*1 

w 

x, Y, Z 

y 

X e , Y e , Z e 

*s> Y s> Z s 


Z 

* 

n 

ft* 


Geodetic latitude of the point S’ 

Successive approximations to 0 
Geocentric latitude of S’ 

Libration phase angle 
Geopotential function 

Coordinates of the observation point P on the geoid surface 

Local (sample) Function 

Earth-fixed coordinates of the satellite 

Tracking station error parameters 

Vector of differential corrections 

Geodetic latitude of point Q; angular arc surface distance 
Geodetic longtitude of point Q 
Astronomical longtitude of P 


oo 


o 


Satellite orbital rate 



REFERENCES 


1. Strange, W. E. , "Marine Geodesy, "EOS, Transactions of the American 
Geophysical Union, March 1971, Vol. 52, No. 3. 

2. Hudson, E. F. , "Determining the Geometric Shape of the Geoid in Ocean- 
Covered Areas of the World by Satellite Radar Altimeter" (abstract), EOS, 
Transactions of the American Geophysical Union, 1970, Vol. 51, p. 263. 

3. New York University School of Engineering and Science, Department of 
Meteorology and Oceanography, Report on Contract N62306-1589, U. S. 

Naval Oceanographic Office, Radar Altimetry from a Spacecraft and its 
Po tential Applications to Geodesy and Oceanography , J. A. Greenwood, 
et al. , May 1967. 

4. Smithsonian Astrophysical Observatory, Report No. 294, Possible Geo- 
potential I mprovement from Satellite Altimetry, C. A. Lundquist, G. E. O. 
Giacaglia, K. Hebb, and S. G. Mair, 1970. 

5. MIT Measurement Systems Lab., RE-37, Combining Satellite Altimetry 
and S urface Gravimetry in Geodetic Determinations , R. G. E. Young, 1970. 

6. Koch, K. R. , "Gravity Anomalies for Ocean Areas from Satellite Altimetry," 
Proceedings of the Second Marine Geodesy Symposium , Marine Technology 
Society, Washington, D. C. , 1970. 

7. National Aeronautics and Space Administration, Goddard Space Flight Center, 
X-552-70-136, Tracking Accuracy Studies of GEOS-C Orbits for Altimetry 
Using Radar and Optical Data , J. G. Marsh and B. C. Douglas, February 
1970 (preprint). 

8. Computer Sciences Corporation, "Task Assignment #210, Determination 

of the Geoid from Satellite Altimeter Data, Memorandum of Understanding," 
R. D. Brown, Letter attachment from R. L. Taylor (CSC) to W. D. Kahn 
(NASA, GSFC), November 4, 1971. 

9. Heiskanen, W. A., and H. Moritz, Physical Geodesy , W. H. Freeman 
and Co. , 1967. 

10. Smithsonian Astrophysical Observatory, Special Report No. 200, Geo- 
detic Parameters for a 1966 Smithsonian Institution Standard Earth , Vol. 

3, C. A. Lundquist and G. Veis, 1966. 


R-l 



11. Richardson, J. A. , " Summary of Activities on Task 256 , Computer Sciences 
Corp. , Technical Memorandum 5035-25600-OITM Task Assignment 256, 
Contract No. NAS5-11790, June 1972. 

12. Vincent, S., W. E. Strange, and J. G. Marsh, "A Detailed Gravimetric 
Geoid from North America to Europe" (paper presented at the National 
Fall Meeting, American Geophysical Union, San Francisco, California, 
December 6-9, 1971). 

13. MIT Experimental Astronomy Laboratory, RE-48, Signal Processing in 
a Satellite Radar Altimeter , C. F. Price, August 1968. 

14. Raytheon Co. , Space Geodesy Altimetry Verification Experiment Design 
Study, Final Report, NASA Contract NASW-1920, E. Weiss, April 1970. 


15. Kershner, R. B., and R. E. Fischell, "Gravity Gradient Stabilization of 
Earth Satellites," Johns Hopkins University Applied Physics Laboratory 
(paper presented at the International Federation of Automatic Control 
Symposium in the Peaceful Uses of Space, Stavanges, Norway, 1965). 

16. Bomford, G. , Geodesy, 2nd Edition. London: Oxford University Press, 
1962. 

17. National Aeronautics and Space Administration, Goddard Space Flight 
Center, Contracts NAS 5-10588 and NAS 5-11729, Mod. 1, Network Ana- 
lysis Program - NAP-3, Mathematical Analysis Documentation , J. J. Lynn, 
J. G. Hartwell, and R. S. Nankervis, June 1, 1970. 

18. National Aeronautics and Space Administration, Goddard Space Flight 
Center, Proposal for Sky Lab. Earth Resources Experimental Package 
Altimeter Evaluation and Analysis , J. W. Bryan, (unpublished). 

19. Johns Hopkins University Applied Physics Laboratory, TG-1024, A Two- 
Quartic Refractivity Profile for the Troposphere, for Correcting Satellite 
Data, H. S. Hopfield, September 1968. 


20. New York University School of Engineering and Science: Department of 
Meteorology and Oceanography, Spacecraft Oceanography Project, U. S. 
Naval Oceanographic Office, Contract N62306-70-A-0075, Task Order 
0001, The Effects of Wind Waves and Swell on the Ranging Accuracy of a 
Radar Altimeter (Part I), W. J. Pierson and E. Mehr. 


R-2 



21. U. S. Department of Commerce, Coast and Geodetic Survey, Special 
Publication No. 98, Manual of Harmonic Analysis and Prediction of Tides , 
P. Shureman, 1958. 

22. National Aeronautics and Space Administration, Goddard Space Flight 
Center, Document 1-500-70-173, The 1975 Satellite Network, September 
1970. 

23. NOAA, Technical Memorandum N059, The Earth's Gravity Field Repre- 
sented by a Simple Layer Potential from Doppler Tracking of Satellites, 

K. R. Koch and B. V. Witte, April 1971. 

24. Vinti, J. P., "Representation of the Earth's Gravitational Potential," 
Celestial Mechanics , December 1971, Vol. 4, No. 3/4. 


R-3 



